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(54) Title: PHASE DETERMINATION OF A RADIATION WAVE FIELD 
(57) Abstract 



A method and apparatus for quantitative de- 
termination of the phase of a radiation wave field 
is disclosed. A representative measure of the rate 
of change of intensity of the radiation wave field 
over a selected surface extending generally across 
the wave field is transformed to produce a first 
integral transform representation. A first filter is 
applied to the first integral transform representa- 
tion corresponding to the inversion of a first dif- 
ferential operator reflected in the measure of rate 
of change of intensity to produce a first modified 
integral transform representation. An inverse of a—* 
the first integral transform is applied to the first 
modified integral transform representation to pro- 
duce an untransformed representation. The un- 
transformed representation is corrected based on lo""*' 
a measure of intensity over said selected surface 
and again transformed to produce a second integral 
transform representation. A second filter is applied 
to the second integral transform representation cor- a— ■** 
responding to the inversion of a second differential 
operator reflected in the corrected untransformed 
representation to produce a second modified inte- 
gral transform representation. An inverse of the 
second integral transform is applied to the second 
modified integral transform representation to pro- 
duce a measure of phase of the radiation wave field 
across the selected plane. 
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PHASE DETERMINATION OF A RADIATION WAVE FIELD 
Field of the Invention 

This invention relates to the determination of phase of a radiation wave field. 
5 The invention also relates to a range of applications in which phase information 
about a radiation wave field can be used. As used in this specification the term 
"radiation wave field" is intended to include all forms of radiation that propagate 
in a wave like manner including but not limited to examples such as X-rays, 
visible light and electrons. 

10 

Background of the Invention 

Techniques for the measurement of the phase of a radiation wave field have 
many applications in fundamental physics and as a basis for a number of 
measurement techniques involving various physical properties. Examples of 
15 applications of phase measurement techniques include the fields of x-ray 
imaging, electron microscopy, optical microscopy as well as optical tomography 
and x-ray phase tomography. 

Phase is typically measured using interferometers of various types. The key 
20 feature of interferometry is the ability to quantitatively measure the phase of a 
wave field. Whilst interferometry based techniques retain significant importance 
it has been recognised that non-interferometric techniques may be used to 
provide phase information. A number of non-interferometric approaches involve 
attempting to solve a transport of intensity equation for a radiation wave field. 
25 This equation relates the irradiance and phase of a paraxial monochromatic 
wave to its longitudinal irradiance derivative and is described in M.R. Teague, 
"Deterministic Phase Retrieval: A Green's Function Solution" J. Opt. Soc. Am. 
73 1434-1441 (1983). The article "Phase imaging by the transport of intensity 
equation" by N. Streibl, Opt. Comm. 49 6-10 (1984), describes an approach 
30 based on the transport of intensity equation by which phase structure can be 
rendered visible by the use of defocus and digital subtraction of intensity data 
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obtained at various defocused distances. This approach only provides for phase 
visualisation and does not provide for the measurement of phase shift. Another 
approach based on solving the transport of intensity equation is disclosed in T.E. 
Gureyev, K.A. Nugent, D. Paganin and A. Roberts, "Rapid phase retrieval using 
5 a Fast Fourier transform", Adaptive Optics, Volume 23, (1995) Optical Society of 
America Technical Digest Series, page 77-79 and T.E. Gureyev and K.A. 
Nugent, "Rapid quantitative phase imaging using the transport of intensity 
equation", Opt. Comm., 133 339-346 (1997). This approach allows the phase of 
a light field to be recovered from two closely spaced intensity measurements 

10 when an illuminating beam has an arbitrary, but everywhere non zero intensity 
distribution limited by rectangular aperture. Whilst this approach can be used for 
non-uniform intensity distributions the extent of the non uniformity is limited and 
introduces significant computational complexity. Consequently this approach is 
not able to cope with non uniformities introduced by some sample absorption 

is profiles or in some intensity illumination distributions. This approach is strictly 
also only applicable to coherent wave fields. 

The article K.A. Nugent, T.E. Gureyev, D.F. Cookson, D. Paganin and Z. Barnea 
"Quantitative phase imaging using hard X-rays" (1996) 77 Phys. Rev. Lett. 2961- 
20 2964 is also based on a solution to the transport of intensity equation. Again the 
technique described cannot be applied to a non-uniform intensity distribution. 

Other approaches based on a solution to the transport of intensity equation 
limited to a requirement of uniformity are described in T.E. Gureyev, K.A. 
25 Nugent, A. Roberts "Phase retrieval with the transport-of-intensity equation: 
matrix solution with the use of Zernike polynomials" J. Opt. Soc. Am. A Vol 12, 
1932-1941 (1995) and T.E. Gureyev, A. Roberts and K.A. Nugent "Partially 
«r coherent fields, the transport-of-intensity equation, and phase uniqueness", J. 

Opt. Soc. Am. A Vol 12, 1942-1946 (1995). 
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A technique for recovery of phase in the case of non-uniform illumination is 
described in T.E. Gureyev and K.A. Nugent "Phase retrieval with the transport- 
of-intensity equation. II. Orthogonal series solution for nonuniform illumination", 
J. Opt. Soc. Am. A Vol 13, 1670-1682 (1996). This approach is based on a 
5 method of orthogonal expansions and can be computationally complex in 
implementation. In many applications this complexity makes the technique 
impractical. 

Disclosure of the Invention 

10 The present invention provides a non-interferometric method and apparatus for 
the measurement of phase. In combination with a direct measurement of 
intensity a measurement of phase allows the phase and intensity at any other 
plane in the radiation wave field to be determined using known techniques. The 
invention also provides the basis for a number of measurement techniques. 

15 

In accordance with a first aspect of the invention there is provided a method of 
quantitative determination of the phase of a radiation wave field including the 
steps of 

(a) producing a representative measure of the rate of change of intensity 
20 of said radiation wave field over a selected surface extending generally across 

the wave field; 

(b) producing a representative measure of intensity of said radiation 
wave field over said selected surface; 

(c) transforming said measure of rate of change of intensity to produce a 
25 first integral transform representation and applying to said first integral transform 

representation a first filter corresponding to the inversion of a first differential 
operator reflected in said measure of rate of change of intensity to produce a 
first modified integral transform representation; 

(d) applying an inverse of said first integral transform to said first 
30 modified integral transform representation to produce an untransformed 

representation; 
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4 

(e) applying a correction based on said measure of intensity over said 
selected surface to said untransformed representation; 

(f) transforming the corrected untransformed representation to produce a 
second integral transform representation and applying to said second integral 
transform representation a second filter corresponding to the inversion of a 
second differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform representation; 

(g) applying an inverse of said second integral transform to said second 
modified integral transform representation to produce a measure of phase of 
said radiation wave field across said selected plane. 

In accordance with a second aspect of the invention there is provided an 
apparatus for quantitative determination of the phase of a radiation wave field 
including 

(a) means to produce a representative measure of the rate of change of 
intensity of said radiation wave field over a selected surface 
extending generally across the direction of propagation; 

(b) means to produce a representative measure of intensity of said 
radiation wave field over said selected surface; 

(c) processing means to sequentially 

(I) transform said measure of rate of change of intensity to 
produce a first integral transform representation; 

(II) apply to said first integral transform representation a first filter 
corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 

(III) apply an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 
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(IV) apply a correction based on said measure of intensity over said 
selected surface to said untransformed representation; 

(V) transform the corrected untransformed representation to 
produce a second integral transform representation; 

(VI) apply to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform 
representation; and 

(VII) apply an inverse of said second integral transform to said 
second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

The selected surface can take any form that extends across the direction of 
propagation of the radiation including planar, part-spherical and part-cylindrical 
surfaces. 

The first and second integral transforms can be of any suitable type and include 
approximations employed for computational convenience, speed or efficiency. 

The first and second integral transforms are preferably produced using a Fourier 
transform. More preferably, the transform is a Fast Fourier transform. The 
method and apparatus of this invention provide for determination of phase of a 
radiation wave field in a manner that is computationally significantly less 
complex than prior art approaches. This results in significantly lower 
computation times. In some examples computation times improved by many 
orders of magnitude have been achieved. 

The first and second differential operators are preferably second order 
differential operators. In the preferred implementation of the method and 
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apparatus the first filter is substantially the same as the second filter. It is further 
preferred that at least one of the first and second filters includes a correction for 
noise in the representative measure of intensity. 

5 In one form of the invention the first filter can comprise selectively suppressing 
first higher frequencies of the first integral transform representation. In this form 
of the invention the second filter can comprise selectively suppressing second 
higher frequencies of said second integral transform representation. 

o The correction based on the measure of intensity over a selected surface can be 
a nil correction where the intensity variations are less than a predetermined 
selected amount. 

Preferably, the measure of the rate of change of intensity and intensity 
15 distribution over the selected surface are produced from measurements of the 
intensity distribution over at least two surfaces extending across the wave field 
and spaced apart along the direction of propagation of the radiation. In another 
form of the invention the representative measure of rate of change intensity in 
the direction of radiation propagation is produced by obtaining a first 
20 representative measurement over a measurement surface extending across the 
direction of propagation for radiation of a first energy and obtaining a second 
representative measurement over said measurement surface for radiation of a 
second different energy. In the case of X-ray radiation, for example, the change 
in radiation energy can be achieved by changing the X-ray target or by suitable 
25 filtering. 

The selected surface for which measurements of intensity and rate of change of 
intensity are produced is preferably located between two of the spaced apart 
surfaces over which intensity distribution is measured. 
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In the preferred form of the invention the selected surface and spaced apart 
surfaces are planar. It is further preferred that the planes are generally 
perpendicular to the average direction of propagation of the radiation. 

The method and apparatus of this invention can be at least partly implemented 
using a suitably programmed computer. In particular the processing means 
preferably comprises a suitably programmed computer and the steps of the 
method are preferably performed using a suitably programmed computer. In 
such forms of the invention intensity input information may take the form of 
digitised images or data containing information from such images. In other 
implementations of the invention a dedicated Fast Fourier transform chip can be 
employed as at least part of the processing means. 

The representative measure of rate of change of intensity is preferably produced 
by subtraction of representative measurements respectively made at locations 
over the spaced apart surfaces. In the preferred form of the invention the 
representative measures of intensity and rate of change of intensity are obtained 
by sampling measurements at selected locations over the surface. Preferably 
the sampling and measurements are made at locations defining a regular array 
over the surface. This can be readily achieved for example by using a CCD 
(charge coupled device) as the detector. 

In the preferred method the direction of propagation of the radiation wave field is 
selected to be the z direction of a Cartesian co-ordinate system and x and y 
components of phase are produced separately. 

In this Cartesian co-ordinate system where the z direction is the direction of 
propagation of the radiation, the preferred filters are of the form 
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(kt+k 2 y )k y 
(k 2 x +k 2 y ) 2 +ak; 



where 

5 k x , k y are the Fourier variables conjugate to x, y and 

a is a constant determined by noise in the intensity measurements and 
is equat to zero for a no noise case. 

The measure of rate of change of intensity is preferably multiplied by the 
10 negative of the average wave number of the radiation before the integral 
transformation into the Fourier domain. 

Th e representative measure of intensity over the spaced apart surfaces can be 
obtained by imaging of that surface through an appropriate imaging system. 
15 That is, the intensity information may be imaged to a detector rather than 
measured at the surface. 



The method of this invention thus provides for the quantitative and decoupled 
determination of phase and intensity of a radiation wave field at any surface 
20 across the direction of propagation of the radiation. From this phase and 
intensity determination it is possible to calculate the phase and intensity at any 
other surface along the direction of propagation. Accordingly, the invention 
provides the basis for a number of measurement techniques. 



25 In a further aspect of the invention there is provided a method of imaging an 
object including the steps of 

(a) exposing the object to a radiation wave field from a source; 

(b) producing a representative measure of the rate of change of 
intensity over a selected surface extending generally across the 
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wave field on the side of the object remote from incident 
radiation; 

(c) producing a representative measure of intensity of said radiation 
wave field over said selected surface; 

(d) transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying to 
said first integral transform representation a first filter 
corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 

(e) applying an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

(f) applying a correction based on said measure of intensity over 
said selected surface to said untransformed representation; 

(g) transforming the corrected untransformed representation to 
produce a second integral transform representation and 
applying to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform 
representation; 

(h) applying an inverse of said second integral transform to said 
second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 



In a still further aspect of the invention there is provided an apparatus for 
imaging an object including 

(a) a source to irradiate the object with a radiation wave field; 



WO 00/26622 



PCT/AU99/00949 



10 



10 



15 



20 



25 



30 



(b) means to produce a representative measure of the rate of 
change of intensity of said radiation wave field over a selected 
surface extending generally across the wave field; 

(c) means to produce a representative measure of intensity of said 
radiation wave field over said selected surface; 

(d) processing means to sequentially 

(I) transform said measure of rate of change of intensity to 
produce a first integral transform representation; 

(II) apply to said first integral transform representation a first 
filter corresponding to inversion of a first differential 
operator reflected in said measure of rate of change of 
intensity to produce a first modified integral transform 
representation; 

(III) apply an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

(IV) apply a correction based on said measure of intensity 
over said selected surface to said untransformed 
representation; 

(V) transform the corrected untransformed representation to 
produce a second integral transform representation; 

(VI) apply to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected 
untransformed representation to produce a second 
modified integral transform representation; and 

(VII) apply an inverse of said second integral transform to said 
second modified integral transform representation to 
produce a measure of phase of said radiation wave field 
across said selected plane. 
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The radiation used to irradiate the object can be a planar wave field or spherical 
wave field or an arbitrary wave field. If it is desired to reproduce the phase in 
the object plane the phase wave field determined by the above method and 
apparatus is back propagated and the wave field used to irradiate is subtracted. 

5 

The method and apparatus of imaging substantially incorporates the 
determination of phase as disclosed in relation to the first and second aspects of 
the invention. The preferred aspects of the invention described in relation to 
those aspects above are also applicable to the method and apparatus of 
10 imaging. 

It is possible in some applications to use a zero object to image plane distance 
corresponding to contact-imaging with zero propagation distance. 



15 



20 



25 



If desired the object can be reconstructed in the object plane by back 
propagating the intensity and quantitative phase information to numerically 
reconstruct an image of the actual object phase and intensity structure. 

In other forms of the method more than two image plane intensity distribution 
measurements can be made to obtain a better estimate of the rate of change of 
intensity or intensity derivative. In this case one or both of the source to object 
or object to image plane distances is changed and another intensity distribution 
measurement is made. The procedure is repeated until the desired number of 
measurements is made. The measurements provide data to which a function 
can be fitted for the determination of rate of change of intensity. 

The method of imaging an object has particular application to point projection 
microscopy using X-rays, visible light or electrons. 



In another aspect this invention provides a method of phase amplitude imaging 
including the steps of 
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irradiating an object with a radiation wave field; 
focussing radiation from the object through an imaging system 
to an imaging surface extending across the wave field 
propagating from the object; 

producing a first representative measure of intensity 
distribution of radiation over said imaging surface at a first 
focus of the imaging system; 

introducing a change in focus of the image on said imaging 
surface through the imaging system; 

producing a second representative measure of intensity 
distribution over said imaging surface; and 
using said first and second representative measures to 
produce a representative measure of intensity and a 
representative measure of rate of change of intensity over a 
selected surface extending across the wave field; 

transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying 
to said first integral transform representation a first filter 
corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 
applying an inverse of said first integral transform to said first 
modified integral transform representation to produce an 

untransformed representation; 

applying a correction based on said measure of intensity over 
said selected surface to said untransformed representation; 
transforming the corrected untransformed representation to 
produce a second integral transform representation and 
applying to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
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representation to produce a second modified integral transform 
representation; 

(k) applying an inverse of said second integral transform to said 
second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

In yet another aspect of this invention there is provided an apparatus for 
phase amplitude imaging of an object including 

a radiation wave field source to irradiate said object; 
an imaging system to focus radiation from said object to an 
imaging surface extending across the wave field propagating from the 
object; 

means to produce a representative measure of radiation 

intensity over said imaging surface; 

said imaging system including selectively operable means to 

adjust said focus of said radiation to said imaging surface to at least a 

first focus and a second focus; 
processing means to: 

(i) produce a representative measure of intensity and a 
representative measure of rate of change of intensity over a 
selected surface extending across the wave field from 
representative measures of radiation intensity over said image 
surface at said first focus and said second focus; 

(ii) transform said measure of rate of change of intensity to produce 
a first integral transform representation; 

(iii) apply to said first integral transform representation a first filter 
corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 
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(iv) apply an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

(v) apply a correction based on said measure of intensity over said 
selected surface to said untransformed representation; 

(vi) transform the corrected untransformed representation to 
produce a second integral transform representation; 

(vii) apply to said second integral transform representation a second 
filter corresponding to the inversion of a second differential 
operator reflected in the corrected untransformed representation 
to produce a second modified integral transform representation; 
and 

(viii) apply an inverse of said second integral transform to said 
second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

Preferably, the numerical aperture of the irradiating wave field is smaller than 
the numerical aperture of the imaging system. 

Preferably, the imaging surface is a detector. The detector is of any suitable 
form, such as for example a CCD camera. 

Preferably the first focus corresponds to an in focus image at the surface and 
the changed focus to a slightly defocussed image. Either negative or positive 
defocus may be used. The defocus is preferably small so that degradation in 
spatial resolution is minimised. In some applications more than two images may 
be obtained to obtain a better estimate of the rate of change of intensity. 

The method and apparatus for phase amplitude imaging substantially 
incorporates the determination of phase as disclosed in relation to the first and 
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second aspects of the invention. The preferred aspects of the invention 
described in relation to those aspects above are also applicable to the method 
and apparatus of imaging. 

5 In a preferred application the method is used for quantitative phase amplitude 
microscopy. In this case the imaging system is a magnification system. 

In the preferred form of the invention the surface is preferably planar. 

10 The invention will now be further described by way of example only, with 
reference to the drawings in which: 

Figure 1 is a schematic illustration of an arrangement for determination of 

phase where an object is illuminated with (a) plane wave radiation and (b) point- 
15 source radiation; 

Figure 2 is a flow chart showing an implementation of the method of phase 

determination in accordance with an embodiment of this invention; 

Figures 3 (a) to (ff) are simulated images illustrating phase determination for 

plane-wave illumination; 
20 Figures 4 (a) to (m) are a series of images illustrating phase determination and 

back propagation to another image plane; 

Figure 5 is a schematic representation of an arrangement for point projection 
microscopy using the method of this invention; 

Figure 6 is a schematic illustration of an arrangement for quantitative phase 
25 amplitude microscopy using the method of this invention; 

Figure 7 is a schematic drawing of an exemplary system for quantitative phase 
amplitude microscopy according to this invention; 

Figure 8(a) to (d) show intensity images and phase images obtained using the 
system shown in Figure 7; 
30 Figur 9 is a graph showing a comparison of measured and expected phase 
profiles for the fibre of Example 3. 



WO 00/26622 



PCT/AU99/00949 



16 

Figur 10 is a schematic drawing of an exemplary system for three dimensional 
optical phase tomography according to this invention. 
Figure 11 is a schematic enlargement of part of Figure 10; 
Figure 12 is a typical tomographic slice through phase image produced in 
5 Example 4; and 

Figure 13 shows a comparison of reconstructed refractive index distribution with 
known refracture index distribution according to Example 4. 

Figures 1(a) and (b) show a schematic arrangement for phase determination in 
10 accordance with this invention where an object is illuminated by plane-wave 
radiation 2 or point source radiation 2 to produce reflected beams 3. 

At each point in space, an optical beam possesses two properties: intensity and 
phase. Intensity is a measure of the amount of energy flowing through each 
15 point, while phase gives a measure of the direction of the energy flow. 

Intensity may be measured directly, for example by recording an image on film. 
Phase is typically measured using interference with a "reference beam". In 
contrast the present method gives a non-interferometric method for measuring 
20 phase. 

Intensity can be measured over two parallel planes A, B extending across the 
direction of propagation of the wave field on the side remote from the incident 
radiation. 

25 

The present invention determines phase by providing a solution to the transport- 
of-intensity equation: 



30 
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where I is the intensity in the plane, the gradient operator in the plane is denoted 
V ± , k is the wave number of the radiation, and dl/dz is the intensity derivative 
or rate of change of intensity. Note that dl/dz is estimated from the difference of 
the measurements in the planes A & B shown in Figure 1, while the intensity I is 
5 given by the average of the measurements. 

In order to obtain a solution to equation 1 the function A is first defined as: 

(2) V ± A = JVJ. 

10 Thus equation (1) becomes: 

(3) V ± .(V x A) = -kd z I. 

Making use of a standard identity V ± • V x = V x 2 , this may be written: 

15 

(4) V x A = -kdJ 

where V x 2 denotes the two-dimensional Laplacian acting over the surface of the 
image. This equation has the following symbolic solution: 

20 

(5) A = -kV x 2 d 2 I . 

If the gradient operator V ± is applied to both sides of this equation, it becomes: 
25 (6) V x A = -kV ± V ± - 2 d 2 I . 

The defining equation (2) for the function A allows (6) to be transformed into: 



(7) JV ± 0 = -kV ± V ± - 2 d z I. 
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Dividing both sides by I then yields: 

(8) V x( z) = -kr l V ± V ± ' 2 dJ . 

5 

Taking the two dimensional divergence V ± • of both sides of (8), and again 
making use of the standard identity V x • V ± = V x 2 , then (8) becomes: 

(9) v x V = -*v x • [r'v^v^cuj. 

10 

This equation has the following symbolic solution: 

(10) 4> = -kv ± (v x •[/" , v x v x - 2 a z /J). 

15 In order to implement a practical solution to equation (10), the following formulae 
are required. A suitably-well-behaved function f(x f y) may be written in the form 
of a two-dimensional Fourier integral: 

(11) Xx 9 y)= J ] f{k x ,k y y' k ^dk x dk r 

-oo -co 

20 

A 

The function f(k x9 k y ) is called the "Fourier transform" of f(x ( y). 
The x derivative of (1 1) yields: 

25 (12) |;/<*. JO = J ] [ *,/( K,k, ) Y k ^dk x dk r 
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Hence the Fourier transform of — f(x,y) is equal to the Fourier transform of 

dx 

f(x,y) multiplied by ik x . Stated differently, — = 'F'^F , where F denotes Fourier 

dx 

transformation and F" 1 denotes inverse Fourier transformation. Similar 

considerations apply to — f(x,y). 

dy 



d 2 d 2 

5 If the Laplacian V 2 = — =- + — =- of (11) is obtained and similar reasoning 

dx dy 



applied, it follows that V" 2 - -F" V 2 F , where k) = k\ + kl . Thus: 



(13) V- 2 =-F-VF, J^/F^F, ±- = *F'%F. 

dx dy 



10 Here, F denotes Fourier transformation, F" 1 denotes inverse Fourier 
transformation, (k x k y ) are the Fourier variables conjugate to (x,y), and 



k r — k x + k 



15 Equations (13) can be used to rewrite equation (10) in the form 



(14) $=<p {x) +<f> {y \< 



= F^k.-^.F/^F- 1 ^^-^^^ 



In practice division by intensity is only performed if that intensity is greater than 
20 a certain threshold value (eg. 0.1% of the maximum value). 



Division by k r does not take place at the point k r = 0 of Fourier space; 
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instead multiplication by zero takes place at this point. This amounts to taking 
the Cauchy principal value of the integral operator V x " 2 . 

In order to quantitatively measure the phase of object it is . necessary to 
5 incorporate some physical constants into the phase recovery algorithm given in 
Equation (14) relating to the experimental setup in use to quantify the variables 
k Xl k y . This can be done by rewriting equation (14) in the following form suitable 
for implementation using a fast Fourier transform: 



10 ^ = < ^+m 



/ _ 2tt 1 / 1 | i j ^ 

9x A&(NAx) 2 i 2 +j 2 /(x,y) i 2 +j 2 X + 



where '-'^"^-'y 



index the frequent components of F(/ + -/_) where the 



intensity derivative dj(x 9 y)\s obtained by subtracting two images /+ and /. 
separated by a distance 8z , i and j are the pixel numbers on the image, and 
15 using the fact that the Fourier space step size is given by 



Ak= 1 



NAx 



where the image is an N x N array of pixels of size Ax . Thus in addition to 
20 measuring the three intensity distributions it is necessary to know the pixel size 
Ax , defocus distance & and wavelength X in order to make a quantitative 
phase measurement. All of these quantities can be readily determined: the pixel 
size can be determined directly for example from the CCD detector geometry (in 
the case of direct imaging), or by existing techniques for calibrating the 
25 transverse distance scales (in the case of an imaging system), the defocus 
distance can be measured directly, and the spectral distribution of the 
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illumination can be determined either by monochromating the incident field or by 
analysing the spectral distribution of the radiation using existing spectroscopic 
methods. 

5 An example of the phase-retrieval method implementing the solution of equation 
(14) can be represented by the flowchart shown in Figure 2. As shown in Figure 
2 the quantitative determination of phase of a radiation wave field proceeds from 
a set of intensity measurements {l„} over the two spaced apart planes A and B. 

A measurement of central intensity l(x,y) in a selected plane parallel to and 
10 midway between the planes A and B is also obtained. The intensity 
measurements are performed over a defined array on each of the two planes A 
and B and the respective values subtracted to produce a measure of the 
intensity derivative. This value is multiplied by the negative of the average wave 
number. The data are split into two component sets and a fast Fourier transform 
15 is performed to produce the respective x and y components in the Fourier 
domain. A filter is then applied to the Fourier domain representations to 
correspond to the inversion of a differential operator reflected in the 
untransformed representation. The differential operator is represented by d~'V* 

for the x component and d~ l V 2 ± for the y component. An inverse Fourier 

20 transform is then performed on each of the x and y components to produce a 
spatial domain value from which the differential operator has been removed. A 
division by the central intensity l(x,y) obtained by averaging the intensity 
measurements over planes A and B is then performed if the intensity level is 
above a selected threshold level. The resultant data is again Fourier 

25 transformed and multiplied by the same filter to again correspond to the 
inversion of a differential operator reflected in the untransformed data. The 
differential operator is again represented by d~ l V 2 ± for the x component and 

d ~y^l± f° r tlie y component . The resultant components are again inverse 
Fourier transformed and summed to provide a retrieved phase measurement. 



30 
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It will be apparent that in general the method according to this invention can 
proceed from any suitable representative determination of intensity derivative or 
rate of change of intensity over a selected surface extending across the 
propagation direction and the intensity over that same surface. As will be 
explained in various examples these data can be obtained in a variety of ways 
and the method implemented to yield phase of the radiation wave field. 

Rewriting equation (14) with: 



£2 y (Ji x7 k y<,ct) — k y k r 



-2 



0(x,y) = 0 M (x,y) + 4> M (x,y), 



gives 



^ x \x,y) = F-'Cl x (k x ,k y ,a)F 



/(x,y) 



1 




(15) \ 



^W) = F- 1 Q,(* x ,*,,a)F 



1 



/(x.y) 



where: 



$(x, y) denotes the recovered phase, 



F denotes Fourier transformation, and F" 1 denotes inverse Fourier 



transformation 



/(x.y) is the intensity distribution over the plane of interest, 
(x,y) are Cartesian coordinates over the plane of interest, 
(k x ,k y ) are the Fourier variables conjugate to (x,y) 
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k =2tt j 



X is the average wavenumber of the radiation, 



X is the average wavelength of the radiation, 
31/ dz is the estimate for the longitudinal intensity derivative, 
a is the regularization parameter used to stabilize the algorithm when 
5 noise is present. 

As given above, the solution to the transport of intensity equation (1) assumes a 
perfect imaging system. That is, there are no "aberrations" present in the optical 
system used to obtain the intensity data which is fed into the algorithm. Of 
10 course, no imaging system is perfect. The imperfections present in an imaging 
system may be quantified by a set of numbers: 



15 which are termed aberration coefficients. 

If intensity data were taken on an imperfect instrument whose imperfections 
were characterized by a certain set of known aberration coefficients A l9 A 29 A 3 ,..., 

it would be desirable if the filters Cl x (k x k y9 d) and Cl y (k x k y9 a) present in (15) 

20 could be replaced by modified filters which explicitly depend upon the aberration 
coefficients: 



25 This would allow the imperfections of the imaging system to be explicitly taken 
into account, leading to quantitatively correct phase retrieval using aberrated 
imaging systems. For the special case of a non-absorbing phase object in a 
radiation wave field of uniform intensity with weak (i.e. much less than 



(16) 



A x , A 2 , A$,... 



(17) 



Cl x (k x k y ,a,A ] ,A 2 ,A 3 ,...) and Cl y {k x k y ,a,A l ,A 2 ,A 3 ,...) 



WO 00/26622 



PCT/AU99/00949 



24 



2n radians) phase variations the appropriate modified filters lead to the following 
functional form for the phase-retrieval algorithm: 

Hp\ a/ v ,A_p-lJ^ F {l aberrated ( X >y)-^\ - [ 

5 where: 

/aberrated(x,y) is the aberrated intensity measured at defocus distance 

&, 

Ann are the aberration coefficients which characterize the imperfect 
imaging system. 



10 



15 



If a filter is defined: 



(19) a.(k x k y ,a J A^A2^^..) = 

1 

^-2x.&J(k 2 x +*, 2 )-2Z m X„ A mn k:K 

Then (18) becomes: 

(20) ^^) = F" 1 QFlF- 1 QF{l aberated (^ } 

The term { I aberrated (x,y)-l }is a measure of rate of change of intensity, lo 

20 intensity is a measurable constant for uniform intensity so that (20) is the same 
general form as (15). Consequently the special case of aberration can be dealt 
with by changing the filter in the general method described above. The x and y 
component filters Q x and Q y are given by 



25 (21) Q x =Q,=-~=rQ 
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Example 1 - Simulations with normally incident plane wave irradiation 

A simulation was conducted in accordance with the arrangement shown in 
5 Figure 1(a) corresponding to planar illumination. The example shows the 
operation of the method on simulated noise-free data. Diffraction patterns are 
calculated using the "angular-spectrum" formalism, an orthodox procedure. 
Figures 3(a) to 3(f) show images produced in the simulation. 

10 Dimensions of all images are 1.00 cm square and provide a sampling array of 
256 x 256 pixels in a plane extending perpendicularly across the propagation 
direction of the radiation. The wavelength of the light was taken to be 632. 8nm. 
The intensity in the plane 2 = 0, which varies from 0 to 1 in arbitrary units, is 
shown in Figure 3(a). Within the area of nonzero illumination, the minimum 

15 intensity was 30% of the maximum intensity. (The black border around the edge 
of the intensity image corresponds to zero intensity.) The input phase, which 
varies from 0 to n radians, is shown in Figure 3(b). 

Images corresponding to planes negatively and positively displaced 2mm from 
20 the z=0 plane shown are in Figures 3(c) and (d) respectively, and have 
respective maximum intensities of 1.60 and 1.75 arbitrary units; the propagation- 
induced phase contrast is clearly visible in each of these images. The two 
defocused images are subtracted to form the intensity derivative, which is shown 
in Figure 3(e). 

25 

Images shown in Figures 3(a) and (e) respectively providing measures of 
intensity and rate of change of intensity across the z=0 plane were then 
processed according to a computer implementation of the method shown in 
Figure 2, in order to yield the recovered phase map given in Figure 3(f). Note 
30 that Figures 3(b) and (f) are plotted on the same greyscale levels, indicating that 
the recovered phase is quantitatively correct. 
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Figures 4 (a) to (h) shows a series of simulated images illustrating phase 
determination and subsequent back-propagation to another image plane. All 
images are 256 pixels x 256 pixels = 1 cm x 1cm in dimensions, with the 

5 radiation wavelength equal to 632.8 nm. The intensity and phase of the 
radiation in a given plane are shown in Figures 4 (a) and (b) respectively. 
Figures 4(c) through (e) respectively show the propagated intensity at 
propagation distances of 199,200 and 201 mm; note the intermixing of 
information from Figures 4 (a) and (b) in the intensity measurements of Figures 

10 4 (c), (d) and (e). Using the images of Figures 4 (c), (d) and (e) only, the phase- 
retrieval algorithm obtained the phase map given in Figure 4(f) for the phase of 
the propagated field at distance 200 mm. Images of Figures (d) and (f) were 
used to numerically back-propagate the field back to the initial plane. This gave 
Figures 4(g) and (h) for the back-propagated intensity and phase, respectively. 

15 These are in excellent agreement with Figures 4 (a) and (b), thus demonstrating 
the use of the phase retrieval techniques for the quantitative determination of the 
amplitude and phase of a field over regions far displaced from those over which 
intensity measurements are made. Note also that the back-propagation is not 
restricted to free space; back-propagation can also be effected through a known 

20 optical system. 

Example 2 - Point projection microscopy 

As shown in Figure 5, radiation such as X-rays, visible light or electrons from a 
25 point source 10 is allowed to propagate through free space to the object 11, 
located at a distance d so from the source. The radiation passes through the 
object 1 1 , and is allowed to propagate a further distance d 0C i to one of the image 
planes h, h. .In in which the intensity of the radiation is detected. This detection 
is performed using a standard device such as a CCD camera, image plate or 
30 other device capable of registering and digitising the intensity distribution. One 
or both of the distances d so and/or d sd is then changed so as to introduce 
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defocus into the images and the intensity distribution is measured once again. 
The case of d O d=0 corresponding to contact-imaging with zero propagation 
distance is included as one possible measurement. 

5 The intensity data is then processed using the above phase recovery method, to 
recover the decoupled intensity and phase information in the imaging plane. 
Parameters, such as wavelength, pixel size, and defocus distances are inserted 
into the algorithm as explained above, to yield quantitative information about the 
magnitude of the phase shift in the image plane. 

io In certain cases a reconstruction of the object in the object plane, as opposed to 
the downstream diffraction planes li...l n , is desired. In this case the intensity and 
quantitative phase information obtained above can be used to back propagate 
the light field to the object plane, thereby numerically reconstructing an image of 
the actual object phase and intensity structure. This can be done using 

15 standard diffraction code. 

In some cases it is desirable to take more than two images in order to obtain a 
better estimate of the intensity derivative dl/dz, in which case one or both of the 
distances d so and/or d sc i is altered once again and another image taken, with 
20 this procedure repeated until the number of desired images is acquired. A 
function can then be fitted to this data from which dl/dz can be computed and 
used in the phase recovery algorithm in place of the simple subtraction of two 
images normally used. 

25 Example 3 - Quantitative phase amplitude microscopy 

Figure 6 schematically shows an arrangement for quantitative phase amplitude 
microscopy. A sample is illuminated using a source of white light Kohler 
illumination 15, commonly found on optical microscopes. The light is transmitted 
30 through an object 16 and collected by a microscope imaging system 17 and 
relayed to a CCD camera 18 or other digital imaging device having a planar 
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imaging surface. Three images are collected: an in-focus image, I 0 , and two 
slightly out of focus images U and I.. The defocus is obtained by suitable means 
such as a drive system 19 to adjust the microscope focus knob. The defocus 
introduced is usually quite small so that degradation in spatial resolution is 
5 minimised, although the optimal amount of defocus to use is determined by 
sample properties and imaging geometry such as magnification, numerical 
apertures, etc. 

When taking the images the numerical aperture of the condenser is chosen to 
10 be less than the numerical aperture of the objective being used. If this is not the 
case then serious image degradation will occur, although the precise amount by 
which the condenser and objective numerical apertures should differ involves a 
tradeoff between image fidelity and spatial resolution, with the optimal difference 
depending on the sample properties and the optics used. 

15 

Intensity data from the collected images 1+ and L are subtracted to produce a 
representative measure of rate of change of intensity (intensity derivative). 
From this and the intensity data of collected image l Q the method described 
above can be used to produce quantitative information about the magnitude of 
20 the phase shift in the image plane. 

As in Example 2 for point projection, there may be cases in which it is desirable 
to take more than two images in order to obtain a better estimate of the intensity 
derivative dl/dz. A function can then be fitted to this data from which dl/dz can 
25 be computed and used in the phase determination method in place of the simple 
subtraction of two images normally used. 

It is also possible to operate this system in a reflection geometry to obtain 
surface topography. The principle of operation is the same, however the optics 
30 have to be folded back on themselves to form a reflection geometry - otherwise 
the process is identical. 
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For certain applications it can also be desirable to filter the light to a particular 
wavelength, although this is not necessary for the described imaging process as 
it works equally well with white light. 

5 

An experimental implementation is shown in Figure 7. An Olympus BX-60 
optical microscope 20 was equipped with a set of UMPIan metallurgical 
objectives and a universal condenser to provide Kohler illumination. In order to 
be able to compare the results with existing imaging modes Nomarski DIC optics 

10 and a set of cover-slip corrected UplanApo objectives were also acquired for this 
microscope, enabling images to be taken of the same field of view using both 
phase retrieval and Nomarski DIC for the purposes of qualitative comparison. A 
12-bit scientific grade Photometries SenSys CCD camera 21 equipped with a 
1300x1035 pixel Kodak KAF-1400 CCD chip was added to the 0.5x video port 

is on the microscope to acquire digital images of the sample. 

The phase recovery technique of this embodiment of the invention requires the 
acquisition of defocused images. A stepper motor drive system 22 was attached 
to the focus knob of the microscope. This stepper motor 22 was attached to the 
20 parallel port of a 133MHz Pentium PC 23 also used to control the CCD camera 
21 , enabling full automation of the acquisition of through-focus image 
sequences. This data acquisition system was linked to custom software written 
to recover phase images from the CCD images, thereby enabling full automation 
of the image acquisition and data processing sequences. 

25 

In order to demonstrate that phase imaging using this invention can accurately 
measure the phase structure of microscopic samples it was necessary to have a 
sample with a well-characterised geometry and refractive index distribution. A 
3M F-SN-3224 optical fibre (a commercial fibre made by 3M) was used. 
30 Independent measurements of the refractive index distribution obtained using 
atomic force microscopy and commercial profiling techniques were available 



WO 00/26622 J PCT/AU99/00949 

30 

enabling accurate prediction of the phase structure of the exit wave field. 
Another advantage of this fibre was that it had three regions of different 
refractive indices, an inner and outer cladding as well as the core, whereas most 
fibres simply have a cladding and core. This provided an additional test for the 
5 phase imaging system because it had to accurately image three transitions in 
refractive index rather than just two. 

The optical fibre was imaged side-on so as to obtain a projection of the 
refractive index through all layers of the fibre structure. This was done by first 

10 stripping the plastic sheath from the fibre by soaking it in isopropyl alcohol then 
using a commercial fibre stripper to remove the plastic coating. A small 
segment of fibre, typically a strand of approximately one to two centimetres in 
length, was placed on a microscope slide, immersed in a pool of index matching 
fluid and covered with a 0.15mm thick cover glass. Any tilt on the cover glass 

is would introduce a spurious tilt into the recovered phase so two small sections of 
fibre, both of similar diameter to the sample, were placed parallel to and about 
0.5cm either side of the main fibre. The cover class was then placed across all 
three fibres to ensure that it was as parallel to the microscope slide as practically 
possible. 

20 

Images of the fibre were taken using an Olympus 40x 0.7NA UplanApo 
objective, which meant that a 500x500 pixel image conveniently spanned the 
whole width of the fibre, and the condenser was set at NA=0.2. The fibre's 
refractive index profiles were known for 632.8nm (HeNe laser) light, so a 

25 625±10nm band-pass interference filter was inserted into the illumination 
system to ensure that the recovered phase profiles were obtained at a 
wavelength as close as possible to that for which data was available on the 
fibre. An intensity image of this sample in the plane of best focus and at ±2//m 
either side of best focus is shown in Figure 8, alongside a phase image 

30 recovered from the two defocused images using the phase-retrieval algorithm 
described above. Note that the fibre is virtually invisible in the in-focus image 
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and barely visible in slightly defocused images, whilst both the fibre and the 
regions of different refractive index, including the 4//m diameter core, are 
clearly visible in the phase image. 

5 Figure 9 shows a comparison of the measured and expected phase profiles with 
the uncertainties indicated in the figure representing one standard deviation of 
the data along the length of the fibre. This variation is thought to be primarily 
due to spatial variations in thickness of the cover glass and microscope slide. 
As can be seen, the recovered and expected phase profiles are in good 
10 agreement with one another, with the predicted profile lying within the error bars 
of the profile produced using the technique of this invention. 

Example 4 - Three-dimensional optical phase tomography 

15 This example demonstrates the application of quantitative phase microscopy to 
the three-dimensional imaging of objects through the use of computed- 
tomography techniques. This is possible using the techniques of this invention 
because the phase shifts introduced by the object can be directly measured 
independently of any intensity variations in the object, thus an inverse Radon 

20 transform can be used to recover the three-dimensional structure directly from 
the projection data. Although the experimental demonstration provided is in the 
optical regime, the same principles are equally applicable to X-ray, electron and 
neutron phase tomography. 

25 For the purposes of collecting three dimensional data sets the same optical 
microscope described in the previous example is used with the addition of a 
rotation stage 24 for the purposes of rotating the sample within the confines of 
the optical microscope imaging area as shown in Figure 10. The rotation stage 
24 is shown in greater detail in Figure 1 1 . 



30 
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The previously described arrangements included a stepper motor drive system 
22 attached to the parallel port of the same 133MHz Pentium PC used to control 
the CCD camera 21 to drive the focus knob of the microscope. A second 
stepper motor 25 was connected to the second channel of the motor drive 
5 system 24 for the purposes of rotating the sample. This data acquisition system 
was linked to custom software written to recover phase images from the CCD 
images, thereby enabling full automation of the image acquisition and data 
processing sequence. Each data set was collected using the same microscope 
as in Example 3 - an Olympus BX-60 optical microscope equipped with a set of 
10 cover-slip corrected UplanApo objectives and a universal condenser to provide 
Kohler illumination. Digital images were captured using a 12-bit Photometries 
SenSys CCD camera equipped with a Kodak KAF-1400 1300x1035 pixel CCD 
chip on the 0.5x video port of the microscope. 

is To prepare the fibre sample 26 for imaging the plastic sheath was removed from 
a small segment of the end of a section of fibre, typically about one centimetre in 
length, by soaking the fibre in isopropyl alcohol then using a commercial fibre 
stripper to remove the plastic coating. The fibre was then cut into a small 
segment of approximately one inch in length, with the unstripped end then being 

20 slid into the end of a 26 gauge, 100mm syringe needle 27 and fixed into position 
with a small amount of 5 minute Araldite. A mount 28 was used to attach the 
needle 27 to stepper motor 25. A pool of index-matching fluid 29 surrounds the 
fibre 26 as shown in Figure 11, with a microscope slide 30 affixed underneath 
the fibre using silicone grease and a 0.15mm thick cover glass 31 placed over 

25 the top. 

Transmission intensity images were collected in the same way as described in 
Example 3 above using an Olympus 20x 0.45NA UMPIan objective with the 
condenser set at NA=0.1. The images taken were 500x500 pixels in size which 
30 conveniently spanned not only the width of the fibre but also the whole region of 
precession. As the refractive index profile of this fibre was known for 632. 8nm 
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(HeNe laser) light, a 625 ± 10nm band-pass interference filter was inserted into 
the illumination system to ensure that the recovered phase profiles were 
obtained at a wavelength as close as possible to that for which data on the fibre 
was available. Each phase image was processed from images taken at ±2/zm 
5 either side of best focus, and data was collected from 100 independent angles 
through 180 degrees equally spaced in steps of 1 .8 degrees between images. A 
typical tomographic phase image is shown in Figure 12. 

The projection data, in the form the reconstructed phase images, were then 

10 processed into three-dimensional data sets using a simple slice-by-slice 
implementation of the summation of filtered back-projections algorithm, with 
code to perform the tomographic reconstruction written in the IDL/PV-Wave 
programming language. First, the data sets were aligned to a common rotation 
axis by taking profiles through the phase data sets and compiling them into a 

is sinogram. A sinusoid was then fitted to prominent features on the data in order 
to determine the location of the rotation axis and the data was digitally shifted so 
that the rotation axis coincided with the middle column of the sinogram to 
simplify the reconstruction process. Fitting a curve to the phase profiles also 
enabled misaligned data sets to be moved back into line, which in turn 

20 improved the quality of the reconstructed image. This realigned projection data 
was then transformed into a single slice through the object by back-projecting 
the collated phase data after filtering the projections to suppress the 1/r point 
spread function associated with back-projected reconstructions. These slices 
through the object were then stacked up on top of each other to produce a full 

25 three-dimensional refractive index distribution of the sample. 

A slice through the reconstructed refractive index distribution is shown in Figure 
13. Note that all three regions of different refractive index are clearly resolved 
and that these regions form concentric cylinders, as is expected for this sample. 
30 A line profile through the centre of this reconstruction is shown in Figure 13 
(dashed line) alongside the known refractive index distribution for this fibre (solid 
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line). The values in the tomographic reconstruction are very close to those of 
* the known profile, which confirms the quantitative phase tomography technique. 
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CLAIMS 

1. A method of quantitative determination of the phase of a radiation 

wave field including the steps of 

(a) producing a representative measure of the rate of change of intensity 
5 of said radiation wave field over a selected surface extending 

generally across the wave field; 

(b) producing a representative measure of intensity of said radiation 
wave field over said selected surface; 

(c) transforming said measure of rate of change of intensity to produce a 
10 first integral transform representation and applying to said first 

integral transform representation a first filter corresponding to the 
inversion of a first differential operator reflected in said measure of 
rate of change of intensity to produce a first modified integral 
transform representation; 

15 (d) applying an inverse of said first integral transform to said first 

modified integral transform representation to produce an 
untransformed representation; 
(e) applying a correction based on said measure of intensity over said 

selected surface to said untransformed representation; 

20 (f) transforming the corrected untransformed representation to produce a 

second integral transform representation and applying to said second 
integral transform representation a second filter corresponding to the 
inversion of a second differential operator reflected in the corrected 
untransformed representation to produce a second modified integral 

25 transform representation; 

(g) applying an inverse of said second integral transform to said second 

modified integral transform representation to produce a measure of 
phase of said radiation wave field across said selected plane. 
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2. A method as claimed in claim 1 wherein said first and second integral 
transforms are produced using a Fourier transform. 

3. A method as claimed in claim 2 wherein said Fourier transform is a 
Fast Fourier transform. 

5 4. A method as claimed in any one of claims 1 to 3 wherein said first 

and second differential operators are second order differential 
operators. 

5. A method as claimed in any one of claims 1 to 4 wherein said first 

filter is substantially the same as said second filter. 
10 6. A method as claimed in any one of claims 1 to 5 wherein said first 

filter includes selectively suppressing first higher frequencies of the 
first integral transform representation. 

7. A method as claimed in any one of claims 1 to 5 wherein at least one 
of said first and second filters includes a correction for noise in said 

15 representative measure of intensity. 

8. A method as claimed in any one of claims 1 to 7 including the step of 
producing said representative measures of intensity and rate of 
change of intensity over said selected surface by producing 
representative measurements corresponding to intensity over at least 

20 two spaced apart surfaces extending across the wave field. 

9. A method as claimed in claim 8 wherein said selected surface is 
between two of said spaced apart surfaces. 

10. A method as claimed in claim 8 wherein said selected surface is one 
of said spaced apart surfaces. 

25 11. A method as claimed in any one of claims 8 to 10 including the step 

of directly detecting representative measures of intensity over said 
spaced apart surfaces. 
12. A method as claimed in any one of claims 8 to 10 including the step 

of producing said representative measure of intensity over at least 

30 one of said spaced apart surfaces by imaging that surface. 
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13. A method as claimed in any one of claims 8 to 12 wherein said 
spaced apart surfaces are substantially parallel. 

14. A method as claimed in claim 13 wherein said spaced apart surfaces 
are substantially planar. 

5 15. A method as claimed in any one of claims 8 to 14 wherein said 

representative measure of rate of change of intensity is produced by 
subtraction of representative measurements of intensity respectively 
made at locations over said spaced apart surfaces. 

16. A method as claimed in any one of claims 1 to 15 wherein said 
io representative measures of intensity and rate of change of intensity 

are obtained by sampling measurements at selected locations over 
said surface. 

17. A method as claimed in claim 16 wherein said sampling 
measurements are made at locations defining a regular array over 

15 said surface. 

18. A method as claimed in claim 2 or claim 3 wherein said radiation 
wave field propagates in a z-direction of a Cartesian co-ordinate 
system and further including the step of producing an x component 
and a y component of phase separately. 

20 19. A method as claimed in claim 18 wherein said first and said second 

filters have a component Q^for producing the x component of phase 
and a component Cl y for producing the y component of phase of the 
form 



25 



(k 2 +k 2 y ) 2 +ak 2 



(k 2 +k 2 y ) 2 +ak 2 y 



where k x k y are the Fourier variables conjugate to x and y; 
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a is a constant determined by noise in the intensity measurements. 

20. A method as claimed in claim 19 including the step of multiplying said 
representative measure of rate of change of intensity by the negative 

5 of average wave number of the radiation before said integral 

transformation. 

21 . A method as claimed in any one of claims 1 to 6 including the step of 
obtaining said representative measure of rate of change of intensity 
by obtaining a first representative measurement over a measurement 

io surface across the wave field for radiation of a first energy and 

obtaining a second representative measurement over said 
measurement surface for radiation of a second different energy. 

22. A method as claimed in claim 2 or claim 3 wherein at least one of said 
first filter and said second filter include a correction for aberrations in 

is said representative measures of intensity and rate of change of 

intensity by including at least one component dependent on the 
aberration, coefficients of a system producing the representative 
measures. 

23. A method as claimed in claim 22 wherein said radiation wave field 
20 propagates in a z-direction of a Cartesian co-ordinate system and 

further including the step of producing an x component and a y 
component of phase separately. 

24. A method as claimed in claim 23 wherein said first and said second 
filters have a component Q^for producing the x component of phase 

25 and a component for producing the y component of phase both of 

the form 

1 

7- An&Mkl + k 2 y ) - 4Z m Z„ A mn Kk n y 
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where k x ,k y are the Fourier variables conjugate to x and y; 
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X is the average wave length of the radiation; 

/aberrated(x,y) is the aberrated intensity measured at defocus distance 
&. 

A mTl are the aberration coefficients which characterize the imperfect 
5 imaging system. 

25. A computer program to execute the steps of any one of claims 1 to 
24. 

26. A computer program stored on computer readable storage media 
including means to execute the steps of any one of claims 1 to 24. 

10 27. An apparatus for quantitative determination of the phase of a 

radiation wave field including 
(a) means to produce a representative measure of the rate of change of 

intensity of said radiation wave field over a selected surface 
extending generally across the wave field; 
15 (b) means to produce a representative measure of intensity of said 

radiation wave field over said selected surface; 
(c) processing means to sequentially 

(I) transform said measure of rate of change of intensity to 
produce a first integral transform representation; 
20 (II) apply to said first integral transform representation a first filter 

corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 

(III) apply an inverse of said first integral transform to said first 
25 modified integral transform representation to produce an 

untransformed representation; 

(IV) apply a correction based on said measure of intensity over said 
selected surface to said untransformed representation; 

(V) transform the corrected untransformed representation to 
30 produce a second integral transform representation; 
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(VI) apply to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform 

5 representation; and 

(VII) apply an inverse of said second integral transform to said 
second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

10 28. An apparatus as claimed in claim 27 wherein said first and second 

integral transforms are produced using a Fourier transform. 

29. An apparatus as claimed in claim 28 wherein said Fourier transform is 
a Fast Fourier transform. 

30. An apparatus as claimed in any one of claims 27 to 29 wherein said 
15 first and second differential operators are second order differential 

operators. 

31. An apparatus as claimed in any one of claims 27 to 30 wherein said 
first filter is substantially the same as said second filter. 

32. An apparatus as claimed in any one of claims 27 to 31 wherein said 
20 first filter includes selectively suppressing first higher frequencies of 

the first integral transform representation. 

33. An apparatus as claimed in any one of claims 27 to 31 wherein at 
least one of said first and second filters includes a correction for noise 
in said representative measure of intensity. 

25 34. An apparatus as claimed in any one of claims 27 to 33 including 

means to produce representative measurements corresponding to 
intensity over at least two spaced apart surfaces extending across the 
wave field. 

35. An apparatus as claimed in claim 34 wherein said selected surface is 

30 between two of said spaced apart surfaces. 
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36. An apparatus as claimed in claim 34 wherein said selected surface is 
one of said spaced apart surfaces. 

37. An apparatus as claimed in any one of claims 34 to 36 including 
detector means positioned to directly detect representative measures 

5 of intensity over said spaced apart surfaces. 

38. An apparatus as claimed in any one of claims 34 to 36 including 
detector means to produce said representative measure of intensity 
over at least one of said spaced apart surfaces and imaging means to 
image that surface onto the detector. 

10 39. An apparatus as claimed in any one of claims 34 to 38 wherein said 

spaced apart surfaces are substantially parallel. 

40. An apparatus as claimed in any one of claims 34 to 38 wherein said 
spaced apart surfaces are substantially planar. 

41. An apparatus as claimed in any one of claims 34 to 40 wherein said 
is means to produce said representative measure of rate of change of 

intensity subtracts representative measurements of intensity 
respectively made at locations over said spaced apart surfaces. 

42. An apparatus as claimed in any one of claims 27 to 41 wherein said 
means to produce a representative measure of intensity and said 

20 means to produce a representative measure of rate of change of 

intensity sample at selected locations over said surface. 

43. An apparatus as claimed in claim 42 wherein said samples are made 
at locations defining a regular array over said surface. 

44. An apparatus as claimed in claim 28 or claim 29 wherein said 
25 radiation wave field propagates in a z-direction of a Cartesian co- 
ordinate system and processing means produces an x component 
and a y component of phase separately. 

45. An apparatus as claimed in claim 44 wherein said processing means 
applies said first and said second filters have a component Q x for 

30 producing the x component of phase and a component Ci y for 

producing the y component of phase of the form 
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(k 2 +k 2 y )k y 



(k 2 +k 2 ) 2 +ak 2 y 



46. 

10 

47. 

15 

48. 

20 

25 49. 



where k x ,k y are the Fourier variables conjugate to x and y; 

a is a constant determent by noise in the intensity measurements. 

An apparatus as claimed in claim 37 wherein said representative 

measure of rate of change of intensity is multiplied by the negative of 

the average wave number of the radiation before said integral 

transformation. 

An apparatus as claimed in any one of claims 27 to 31 wherein said 
representative measure of rate of change of intensity is produced by 
obtaining a first representative measurement over a measurement 
surface across the wave field for radiation of a first energy and 
obtaining a second representative measurement over said 
measurement surface for radiation of a second different energy. 
An apparatus as claimed in claim 28 or claim 29 wherein at least one 
of said first filter and said second filter include a correction for 
aberrations in said representative measures of intensity and rate of 
change of intensity by including at least one component dependent 
on the aberration, coefficients of a system producing the 
representative measures. 

An apparatus as claimed in claim 48 wherein said radiation wave field 
propagates in a z-direction of a Cartesian co-ordinate system and 
wherein an x component and a y component of phase are produced 
separately. 
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50. An apparatus as claimed in claim 49 wherein said first and said 

second filters have a component Q^for producing the x component of 

phase and a component Cl y for producing the y component of phase 

both of the form 

5 

1 



where k Xl k y are the Fourier variables conjugate to x and y; 
J is the average wave length of the radiation; 
10 /aberrated(x l y) is the aberrated intensity measured at defocus distance 

&. 

v4 mn are the aberration coefficients which characterize the imperfect 
imaging system. 

51 . A method of imaging an object including the steps of 

is (a) exposing the object to a radiation wave field from a source; 

(b) producing a representative measure of the rate of change of 
intensity over a selected surface extending generally across the 
wave field on the side of the object remote from incident 
radiation; 

20 (c) producing a representative measure of intensity of said radiation 

wave field over said selected surface; 

(d) transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying to 
said first integral transform representation a first filter 

25 corresponding to the inversion of a first differential operator 

reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 

(e) applying an inverse of said first integral transform to said first 
modified integral transform representation to produce an 

30 untransformed representation; 



WO 00/26622 



PCT/AU99/00949 



44 



(f) applying a correction based on said measure of intensity over 
said selected surface to said untransformed representation; 

(g) transforming the corrected untransformed representation to 
produce a second integral transform representation and 

5 applying to said second integral transform representation a 

second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform 
representation; 

10 (h) applying an inverse of said second integral transform to said 

second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

52. A method as claimed in claim 51 including the step of producing said 
15 representative measures of intensity and rate of change of intensity 

over said selected surface by producing representative 
measurements corresponding to intensity over at least two spaced 
apart surfaces extending across the wave field. 

53. A method as claimed in claim 52 wherein said selected surface is 
20 between two of said spaced apart surfaces. 

54. A method as claimed in claim 52 wherein said selected surface is one 
of said spaced apart surfaces. 

55. A method as claimed in any one of claims 52 to 54 wherein said 
spaced apart surfaces are substantially parallel. 

25 56. A method as claimed in any one of claims 52 to 54 wherein said 

spaced apart surfaces are substantially planar. 
57. A method as claimed in any one of claims 52 to 56 wherein said 

representative measure of rate of change of intensity is produced by 
subtraction of representative measurements of intensity respectively 

30 made at locations over said spaced apart surfaces. 



WO 00/26622 



PCT/AU99/00949 



45 



58. A method as claimed in claim 51 including the step of producing said 
measures of intensity and rate of change of intensity over said 
selected surface by producing first representative measurements 
corresponding to intensity over a first surface extending across the 

5 wave field, changing the distance between said source and said 

object, and producing second representative measurements 
corresponding to intensity over said first surface for the changed 
distance between said object and said source. 

59. A method as claimed in claim 58 wherein said selected surfaces is 
io one of said spaced apart surfaces. 

60. A method as claimed in any one of claims 51 to 57 including the step 
of directly detecting representative measures of intensity over said 
spaced apart surfaces. 

61. A method as claimed in any one of claims 51 to 60 wherein said 
is selected surface is spaced apart from said object in the direction of 

propagation of said radiation. 

62. A method as claimed in any one of claims 51. to 61 wherein said 
source is substantially a point source. 

63. A method as claimed in any one of claims 51 to 62 wherein said first 
20 and second integral transforms are produced using a Fourier 

transform. 

64. A method as claimed in any one of claims 51 to 63 wherein said 
Fourier transform is a Fast Fourier transform. 

65. An apparatus for imaging an object including : 

25 (a) a source to irradiate the object with a radiation wave field; 

(b) means to produce a representative measure of the rate of 
change of intensity of said radiation wave field over a selected 
surface extending generally across the wave field. 

(c) means to produce a representative measure of intensity of said 
30 radiation wave field over said selected surface; 

(d) processing means to sequentially 



WO 00/26622 



PCT/AU99/00949 



(i) transform said measure of rate of change of intensity to 
produce a first integral transform representation; 

(ii) apply to said first integral transform representation a first 
filter corresponding to the inversion of a first differential 

5 operator reflected in said measure of rate of change of 

intensity to produce a first modified integral transform 
representation; 

(iii) apply an inverse of said first integral transform to said first 
modified integral transform representation to produce an 

to untransformed representation; 

(iv) apply a correction based on said measure of intensity over 
said selected surface to said untransformed 
representation; 

(v) transform the corrected untransformed representation to 
15 produce a second integral transform representation; 

(vi) apply to said second integral transform representation a 
second filter corresponding to the inversion of a second 
differential operator reflected in the corrected 
untransformed representation to produce a second 

20 modified integral transform representation; and 

(vii) apply an inverse of said second integral transform to said 
second modified integral transform representation to 
produce a measure of phase of said radiation wave field 
across said selected plane. 

25 66. An apparatus as claimed in claim 65 including means to produce 

representative measurements corresponding to intensity over at least 
two spaced apart surfaces extending across the wave field. 
67. An apparatus as claimed in claim 66 wherein said selected surface is 

between two of said spaced apart surfaces. 

30 68. An apparatus as claimed in claim 66 wherein said selected surface is 

one of said spaced apart surfaces. 



WO 00/26622 



PCT/AU99/00949 



47 

69. An apparatus as claimed in any one of claims 66 to 68 including 
detector means positioned to directly detect representative measures 
of intensity over said spaced apart surfaces. 

70. An apparatus as claimed in any one of claims 66 to 68 including 
detector means to produce said representative measure of intensity 
over at least one of said spaced apart surfaces and imaging means to 
image that surface onto the detector. 

71. An apparatus as claimed in any one of claims 66 to 70 wherein said 
spaced apart surfaces are substantially parallel. 

72. An apparatus as claimed in any one of claims 66 to 71 wherein said 
representative measure of rate of change of intensity is produced by 
subtraction of representative measurements of intensity respectively 
made at locations over said spaced apart surfaces. 

73. An apparatus as claimed in claim 65 including means to produce said 
measures of intensity and rate of change of intensity overs said 
selected surface by producing first representative measurements 
corresponding to intensity over a first surface extending across the 
wave field; means to change the distance between said source and 
said object, and means to produce second representative 
measurements corresponding to intensity over said first surface for 
the changed distance between said object and said source. 

74. An apparatus as claimed in claim 73 wherein said selected surface is 
one of said spaced apart surfaces. 

75. An apparatus as claimed in any one of claims 65 to 72 including 
means to directly detecting representative measures of intensity over 
said spaced apart surfaces. 

76. An apparatus as claimed in any one of claims 65 to 75 wherein said 
selected surface is spaced apart from said object in the direction of 
propagation of said radiation. 

77. An apparatus as claimed in any one of claims 65 to 76 wherein said 
source is substantially a point source. 
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78. An apparatus as claimed in any one of claims 65 to 77 wherein said 
first and second integral transforms are produced using a Fourier 
transform. 

79. An apparatus as claimed in claim 78 wherein said Fourier transform is 
5 a Fast Fourier transform. 

80. A method of phase amplitude imaging including the steps of 

(a) irradiating an object with a radiation wave field; 

(b) focussing radiation from the object through an imaging system 
to an imaging surface extending across the wave field 

io propagating from the object; 

(c) producing a first representative measure of intensity 
distribution of radiation over said imaging surface at a first 
focus of the imaging system; 

(d) introducing a change in focus of the image on said imaging 
is surface through the imaging system; 

(e) producing a second representative of measure intensity 
distribution over said imaging surface; and 

(f) using said first and second representative measures to 
produce a representative measure of intensity and a 

20 representative measure of rate of change of intensity in the 

direction of radiation propagation over a selected surface 
extending across the wave field; 

(g) transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying 

25 to said first integral transform representation a first filter 

corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 
produce a first modified integral transform representation; 

(h) applying an inverse of said first integral transform to said first 
30 modified integral transform representation to produce an 

untransformed representation; 
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(i) applying a correction based on said measure of intensity over 

said selected surface to said untransformed representation; 
(j) transforming the corrected untransformed representation to 
produce to said second integral transform representation a 
5 second integral transform representation and applying a 

second filter corresponding to the inversion of a second 
differential operator reflected in the corrected untransformed 
representation to produce a second modified integral transform 
representation; 

io (k) applying an inverse of said second integral transform to said 

second modified integral transform representation to produce a 
measure of phase of said radiation wave field across said 
selected plane. 

81. A method as claimed in claim 80 wherein said radiation wave field 
15 has a numerical aperture smaller than the numerical aperture of said 

imaging system. 

82. A method as claimed in claim 80 or claim 81 said first focus of the 
imaging system produces an infocus image at the imaging surface 
and said second focus of the imaging system produces a slightly 

20 defocused image at the imaging surface. 

83. A method as claimed in any one of claims 80 to 82 wherein said 
imaging surface is substantially planar. 

84. A method as claimed in any one of claims 80 to 83 wherein the 
imaging surface is an intensity detector. 

25 85. A method as claimed in any one of claims 80 to 84 wherein said 

imaging surface is said selected surface. 

86. A method as claimed in any one of claims 80 to 85 wherein said 
integral transform is a Fourier transform. 

87. A method as claimed in claim 86 wherein said Fourier transform is a 
30 Fast Fourier transform. 
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88. A method as claimed in any one of claims 80 to 87 wherein said 

representative measure of rate of change of intensity is produced by 
subtraction of said first and second representative measurements of 
intensity. 

5 89. A method as claimed in any one of claims 80 to 88 wherein said 

representative measures of intensity and rate of change of intensity 
are obtained by sampling measurements at selected location over 
said imaging surface. 

90. A method as claimed in claim 89 wherein said sampling 
I0 measurements are made at locations defining a regular array over 

said imaging surface. 

91. A method as claimed in claim 87 wherein said radiation wave field 
propagates in a z-direction of a cartesian co-ordinate system and 
further including the step of producing an x component and a y 

15 component of phase separately. 

92. A method as claimed in claim 91 wherein said first and said second 
filters have a component Q. x for producing the x component of phase 
and a component O y for producing the y component of phase of the 
form 

<*, 2 



20 Q x = 



{k 2 x +k 2 ) 2 +ak] 



(K+k 2 y )k y 

(k 2 x +k 2 y ) 2 +akl 



where k x ,k y are the Fourier variables conjugate to x and y; 
25 a is a constant determined by noise in the intensity measurements. 

93. A method as claimed in claim 92 including the step of multiplying said 

representative measure of rate of change of intensity by the negative 
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of the average wave number of the radiation before said integral 
transformation. 

94. An apparatus for phase amplitude imaging of an object including 

a radiation wave field source to irradiate said object; 
5 an imaging system to focus radiation from said object to an 

imaging surface extending across the wave field propagating from the 
object; 

means to produce a representative measure of radiation 
intensity over said imaging surface; 
10 said imaging system including selectively operable means to 

adjust said focus of said radiation to said imaging surface to at least 
at a first focus and a second focus; 

processing means to: 
(i) produce a representative measure of intensity and a 
15 representative measure of rate of change of intensity in the 

direction of radiation propagation over a selected surface 
extending across the wave field from representative 
measures of radiation intensity over said image surface at 
said first focus and said second focus; 
20 (ii) transform said measure of rate of change of intensity to 

produce a first integral transform representation; 

(iii) apply to said first integral transform representation a first filter 
corresponding to the inversion of a first differential operator 
reflected in said measure of rate of change of intensity to 

25 produce a first modified integral transform representation; 

(iv) apply an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

(v) apply a correction based on said measure of intensity over 
30 said selected surface to said untransformed representation; 
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(vi) transform the corrected untransformed representation to 
produce a second integral transform representation; 

(vii) apply to said second integral transform representation a 
second filter corresponding to the inversion of a second 

5 differential operator reflected in the corrected untransformed 

representation to produce a second modified integral 
transform representation; and 

(viii) apply an inverse of said second integral transform to said 
second modified integral transform representation to produce 

io a measure of phase of said radiation wave field across said 

selected plane. 

95. An apparatus as claimed in claim 94 wherein said radiation wave field 

has a numerical aperture smaller than the numerical aperture of said 
imaging system. 

15 96. An apparatus as claimed in claim 94 or claim 95 said first focus of the 

imaging system produces an infocus image at the imaging surface 
and said second focus of the imaging system produces a slightly 
defocused image at the imaging surface. 

97. An apparatus as claimed in any one of claims 94 to 96 wherein said 
20 imaging surface is substantially planar. 

98. An apparatus as claimed in any one of claims 94 to 97 wherein the 
imaging surface is an intensity detector. 

99. An apparatus as claimed in any one of claims 94 to 98 wherein said 
imaging surface is said selected surface. 

25 100. An apparatus as claimed in any one of claims 94 to 99 wherein said 
integral transform is a Fourier transform. 

101. An apparatus as claimed in claim 100 wherein said Fourier transform 
is a Fast Fourier transform. 

102. An apparatus as claimed in any one of claims 94 to 101 wherein said 
30 representative measure of rate of change of intensity is produced by 
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subtraction of said first and second representative measurements of 
intensity. 

103. An apparatus as claimed in any one of claims 94 to 102 wherein said 
representative measures of intensity and rate of change of intensity 

5 are obtained by sampling measurements at selected location over 

said imaging surface. 

104. An apparatus as claimed in claim 103 wherein said sampling 
measurements are made at locations defining a regular array over 
said imaging surface. 

10 105. An apparatus as claimed in claim 101 wherein said radiation wave 
field propagates in a z-direction of a cartesian co-ordinate system and 
further including the step of producing an x component and a y 
component of phase separately. 
106. An apparatus as claimed in claim 105 wherein said first and said 

15 second filters have a component fi^for producing the x component 

of phase and a component Q y for producing the y component of 
phase of the form 



(k 2 +k 2 ) 2 +ak 2 



20 



(k 2 +k 2 ) 2 +ak 2 



where k Xl k y are the Fourier variables conjugate to x and y; 
a is a constant determent by noise in the intensity measurements. 
25 107. A method as claimed in claim 106 including the step of multiplying 
said representative measure of rate of change of intensity by the 
negative of the average wave number of the radiation before said 
integral transformation. 
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108. A computer program for processing a representative measure of the 
rate of change of intensity of a radiation wave field over a selected 
surface extending generally across the wave field and a 
representative measure of intensity of said radiation wave field over 
5 said selected surface, said program including: 

(a) code for transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying to said 
first integral transform representation a first filter corresponding to the 
inversion of a first differential operator reflected in said measure of 

10 rate of change of intensity to produce a first modified integral 

transform representation; 

(b) code for applying an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

15 (c) code for applying a correction based on said measure of intensity 

over said selected surface to said untransformed representation; 

(d) code for transforming the corrected untransformed representation to 
produce a second integral transform representation and applying to 
said second integral transform representation a second filter 

20 corresponding to the inversion of a second differential operator 

reflected in the corrected untransformed representation to produce a 
second modified integral transform representation; 

(e) code for applying an inverse of said second integral transform to said 
second modified integral transform representation. 

25 109. A computer program stored on computer readable storage media for 
processing a representative measure of the rate of change of 
intensity of a radiation wave field over a selected surface extending 
generally across the wave field and a representative measure of 
intensity of said radiation wave field over said selected surface, said 

30 program including: 
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(a) code for transforming said measure of rate of change of intensity to 
produce a first integral transform representation and applying to said 
first integral transform representation a first filter corresponding to the 
inversion of a first differential operator reflected in said measure of 
rate of change of intensity to produce a first modified integral 
transform representation; 

(b) code for applying an inverse of said first integral transform to said first 
modified integral transform representation to produce an 
untransformed representation; 

(c) code for applying a correction based on said measure of intensity 
over said selected surface to said untransformed representation; 

(d) code for transforming the corrected untransformed representation to 
produce a second integral transform representation and applying to 
said second integral transform representation a second filter 
corresponding to the inversion of a second differential operator 
reflected in the corrected untransformed representation to produce a 
second modified integral transform representation; 

(e) code for applying an inverse of said second integral transform to said 
second modified integral transform representation. 
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(e) intensity derivative 
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Figure 3 




(a) intensity, z = 0 mm (b) phase, z = 0 mm 




(g) backpropagated intensity, z = 0 mm (h) backpropagated phase, z = 0 mm 

Figure 4 
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